Last updated: 2024-08-20
Checks: 7 0
Knit directory: DCM_snRNAseq/
This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20240606) was run prior to running
the code in the R Markdown file. Setting a seed ensures that any results
that rely on randomness, e.g. subsampling or permutations, are
reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version c979c9e. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for
the analysis have been committed to Git prior to generating the results
(you can use wflow_publish or
wflow_git_commit). workflowr only checks the R Markdown
file, but you know if there are other scripts or data files that it
depends on. Below is the status of the Git repository when the results
were generated:
Ignored files:
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: analysis/figure/
Untracked files:
Untracked: analysis/omnipathr-log/
Untracked: code/Add_metadata.R
Untracked: code/Check DYSF.R
Untracked: code/Clustering_genes.R
Untracked: code/DE_5_percent.Rmd
Untracked: code/DE_5_percent.html
Untracked: code/DE_5_percent/
Untracked: code/DE_CM_test1.R
Untracked: code/DE_no_401.Rmd
Untracked: code/DE_no_401.html
Untracked: code/DE_no_401/
Untracked: code/Differential abundance test1.R
Untracked: code/Differential_Expression_edgeR_All.Rmd
Untracked: code/Differential_Expression_edgeR_All_2.Rmd
Untracked: code/Differential_Expression_edgeR_All_2.html
Untracked: code/Differential_Expression_edgeR_All_Age.Rmd
Untracked: code/Differential_Expression_edgeR_All_Age_2.Rmd
Untracked: code/Differential_Expression_edgeR_All_Age_2.html
Untracked: code/Differential_Expression_edgeR_All_groups.Rmd
Untracked: code/Differential_Expression_edgeR_All_groups.html
Untracked: code/Differential_Expression_edgeR_All_groups_2.Rmd
Untracked: code/Differential_Expression_edgeR_All_groups_2.html
Untracked: code/Differential_Expression_edgeR_All_groups_3.Rmd
Untracked: code/Differential_Expression_edgeR_All_groups_3.html
Untracked: code/PCA and CCA analysis test.R
Untracked: code/Pseudobulk DCM HC.R
Untracked: code/QC_integration_annotation_orig.Rmd
Untracked: code/UpSet_plot_DEGs.Rmd
Untracked: code/Volcano_highlighted_genes.R
Untracked: code/ezInteractiveTable.Rmd
Untracked: code/ezInteractiveTable.html
Untracked: code/old.R
Untracked: core
Untracked: data/Cellbender_output/
Untracked: data/Cellranger_output/
Untracked: data/DCM_Clinical_data.xlsx
Untracked: data/DCM_Clinical_data_26.xlsx
Untracked: data/Raw/
Untracked: omnipathr-log/
Untracked: output/Cardiomyocytes_DA_DE/
Untracked: output/Cardiomyocytes_DA_DE_26/
Untracked: output/Cardiomyocytes_subclustering/
Untracked: output/Cardiomyocytes_subclustering_26/
Untracked: output/Differential_expression_edgeR_All/
Untracked: output/Differential_expression_edgeR_All_Age/
Untracked: output/QC_integration_annotation/
Untracked: output/QC_integration_annotation_26/
Unstaged changes:
Modified: analysis/Cardiomyocytes_DA_DE.Rmd
Deleted: analysis/Differential_Expression_edgeR_All.Rmd
Deleted: analysis/Differential_Expression_edgeR_All_Age.Rmd
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were
made to the R Markdown
(analysis/Cardiomyocytes_DA_DE_26.Rmd) and HTML
(docs/Cardiomyocytes_DA_DE_26.html) files. If you’ve
configured a remote Git repository (see ?wflow_git_remote),
click on the hyperlinks in the table below to view the files as they
were in that past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| Rmd | c979c9e | GinoBonazza | 2024-08-20 | wflow_publish("analysis/Cardiomyocytes_DA_DE_26.Rmd") |
| html | 182f9de | GinoBonazza | 2024-08-19 | Build site. |
| html | 5cd0f74 | GinoBonazza | 2024-08-01 | Build site. |
| Rmd | f3f2075 | GinoBonazza | 2024-08-01 | wflow_publish("analysis/Cardiomyocytes_DA_DE_26.Rmd") |
| html | 61e8d17 | GinoBonazza | 2024-08-01 | Build site. |
| Rmd | 9e4b482 | GinoBonazza | 2024-08-01 | wflow_publish("analysis/Cardiomyocytes_DA_DE_26.Rmd") |
# Get current file name to make folder
current_file <- "Cardiomyocytes_DA_DE_26"
# Load libraries
library(here)
library(readr)
library(readxl)
library(xlsx)
library(Seurat)
library(DropletUtils)
library(Matrix)
library(scDblFinder)
library(scCustomize)
library(dplyr)
library(ggplot2)
library(magrittr)
library(tidyverse)
library(reshape2)
library(S4Vectors)
library(SingleCellExperiment)
library(pheatmap)
library(png)
library(gridExtra)
library(knitr)
library(scales)
library(RColorBrewer)
library(Matrix.utils)
library(tibble)
library(ggplot2)
library(scater)
library(patchwork)
library(statmod)
library(ArchR)
library(clustree)
library(harmony)
library(gprofiler2)
library(clusterProfiler)
library(org.Hs.eg.db)
library(AnnotationHub)
library(ReactomePA)
library(statmod)
library(edgeR)
library(speckle)
library(EnhancedVolcano)
library(decoupleR)
library(OmnipathR)
library(dorothea)
library(enrichplot)
library(png)
library(ezRun)
library(UpSetR)
library(ComplexHeatmap)
#Output paths
output_dir_data <- here::here("output", current_file)
if (!dir.exists(output_dir_data)) dir.create(output_dir_data)
if (!dir.exists(here::here("docs", "figure"))) dir.create(here::here("docs", "figure"))
output_dir_figs <- here::here("docs", "figure", paste0(current_file, ".Rmd"))
if (!dir.exists(output_dir_figs)) dir.create(output_dir_figs)
Load cardiomyocytes dataset
CM <- readRDS(here::here("output", "Cardiomyocytes_subclustering_26", "CM.rds"))
DefaultAssay(CM) <- "RNA"
CM <- NormalizeData(CM)
#Prepare matadata for differential abundance and differential expression testing
metadata <- CM@meta.data %>%
dplyr::select(Sample:CO_l.min) %>%
unique() %>%
dplyr::select(-LVID_mm, -RVIT_mm)
metadata <- metadata %>%
dplyr::select(1:3, 5, 6, 4, 7:ncol(metadata)) %>%
dplyr::arrange(match(Sample, names(table(CM$Sample))))
rownames(metadata) <- metadata$Sample
# Check if metadata$Sample matches names(table(CM$Sample))
all.equal(metadata$Sample, names(table(CM$Sample)))
[1] TRUE
ezInteractiveTableRmd(metadata)
# Filter for Group == "DCM"
metadata_DCM <- metadata %>%
dplyr::filter(Group == "DCM")
# Initialize an empty data frame for storing differential abundance results
differential_abundance <- data.frame()
props <- getTransformedProps(CM$cell_state, CM$Sample, transform="logit")
# Loop through each metadata column from Age_years to CO_l.min
for (i in 6:(length(metadata_DCM))) {
metadata_subset <- metadata_DCM[complete.cases(metadata_DCM[[i]]),]
props_subset <- props
props_subset$Counts <- props$Counts[, metadata_subset$Sample, drop = FALSE]
props_subset$TransformedProps <- props$TransformedProps[, metadata_subset$Sample, drop = FALSE]
props_subset$Proportions <- props$Proportions[, metadata_subset$Sample, drop = FALSE]
metadata_subset <- metadata_subset %>%
arrange(Sample, colnames(props_subset$Proportions))
design <- model.matrix(~ metadata_subset[[i]])
fit <- lmFit(props_subset$TransformedProps, design)
fit <- eBayes(fit, robust=TRUE)
differential_abundance_temp <- topTable(fit, n = Inf, coef = 2)
differential_abundance_temp$metadata <- colnames(metadata_subset[i])
differential_abundance_temp$cell_state <- rownames(differential_abundance_temp)
differential_abundance <- rbind(differential_abundance, differential_abundance_temp)
}
# Clean up the environment
rm(metadata_subset, props_subset, design, fit, differential_abundance_temp)
# Reset row names and arrange the results by adjusted p-value
rownames(differential_abundance) <- NULL
differential_abundance <- dplyr::arrange(differential_abundance, adj.P.Val)
write.csv(differential_abundance, file = here::here(output_dir_data, "CM_Differential_Abundance.csv"), quote=F, row.names = F)
t <- ezInteractiveTableRmd(differential_abundance)
t[["x"]][["options"]][["pageLength"]] <- 10
t
metadata_subset <- metadata[complete.cases(metadata[["LVEDD_mm"]]),]
props_subset <- props
props_subset$Counts <- props$Counts[, metadata_subset$Sample, drop = FALSE]
props_subset$TransformedProps <- props$TransformedProps[, metadata_subset$Sample, drop = FALSE]
props_subset$Proportions <- props$Proportions[, metadata_subset$Sample, drop = FALSE]
metadata_subset <- metadata_subset %>%
arrange(Sample, colnames(props_subset$Proportions))
design <- model.matrix(~ metadata_subset$LVEDD_mm)
fit.prop <- lmFit(props_subset$Proportions,design)
par(mfrow=c(1,5), mar = c(5, 5, 3, 1))
for(i in seq(1,5,1)){
plot(metadata_subset$LVEDD_mm, props_subset$Proportions[i,], main = rownames(props_subset$Proportions)[i],
pch=16, cex=2, ylab="Proportions", xlab="LVEDD_mm", cex.lab=2, cex.axis=1.5,
cex.main=2.5)
abline(a=fit.prop$coefficients[i,1], b=fit.prop$coefficients[i,2], col=4,
lwd=1)
}

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
rm(metadata_subset, props_subset, design, fit.prop)
metadata_subset <- metadata[complete.cases(metadata[["LVEF_percent"]]),]
props_subset <- props
props_subset$Counts <- props$Counts[, metadata_subset$Sample, drop = FALSE]
props_subset$TransformedProps <- props$TransformedProps[, metadata_subset$Sample, drop = FALSE]
props_subset$Proportions <- props$Proportions[, metadata_subset$Sample, drop = FALSE]
metadata_subset <- metadata_subset %>%
arrange(Sample, colnames(props_subset$Proportions))
design <- model.matrix(~ metadata_subset$LVEF_percent)
fit.prop <- lmFit(props_subset$Proportions,design)
par(mfrow=c(1,5), mar = c(5, 5, 3, 1))
for(i in seq(1,5,1)){
plot(metadata_subset$LVEF_percent, props_subset$Proportions[i,], main = rownames(props_subset$Proportions)[i],
pch=16, cex=2, ylab="Proportions", xlab="LVEF_percent", cex.lab=2, cex.axis=1.5,
cex.main=2.5)
abline(a=fit.prop$coefficients[i,1], b=fit.prop$coefficients[i,2], col=4,
lwd=1)
}

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
rm(metadata_subset, props_subset, design, fit.prop)
metadata_subset <- metadata[complete.cases(metadata[["Age_years"]]),]
props_subset <- props
props_subset$Counts <- props$Counts[, metadata_subset$Sample, drop = FALSE]
props_subset$TransformedProps <- props$TransformedProps[, metadata_subset$Sample, drop = FALSE]
props_subset$Proportions <- props$Proportions[, metadata_subset$Sample, drop = FALSE]
metadata_subset <- metadata_subset %>%
arrange(Sample, colnames(props_subset$Proportions))
design <- model.matrix(~ metadata_subset$Age_years)
fit.prop <- lmFit(props_subset$Proportions,design)
par(mfrow=c(1,5), mar = c(5, 5, 3, 1))
for(i in seq(1,5,1)){
plot(metadata_subset$Age_years, props_subset$Proportions[i,], main = rownames(props_subset$Proportions)[i],
pch=16, cex=2, ylab="Proportions", xlab="Age_years", cex.lab=2, cex.axis=1.5,
cex.main=2.5)
abline(a=fit.prop$coefficients[i,1], b=fit.prop$coefficients[i,2], col=4,
lwd=1)
}

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
rm(metadata_subset, props_subset, design, fit.prop)
metadata_subset <- metadata[complete.cases(metadata[["RVSP_mmHg"]]),]
props_subset <- props
props_subset$Counts <- props$Counts[, metadata_subset$Sample, drop = FALSE]
props_subset$TransformedProps <- props$TransformedProps[, metadata_subset$Sample, drop = FALSE]
props_subset$Proportions <- props$Proportions[, metadata_subset$Sample, drop = FALSE]
metadata_subset <- metadata_subset %>%
arrange(Sample, colnames(props_subset$Proportions))
design <- model.matrix(~ metadata_subset$RVSP_mmHg)
fit.prop <- lmFit(props_subset$Proportions,design)
par(mfrow=c(1,5), mar = c(5, 5, 3, 1))
for(i in seq(1,5,1)){
plot(metadata_subset$RVSP_mmHg, props_subset$Proportions[i,], main = rownames(props_subset$Proportions)[i],
pch=16, cex=2, ylab="Proportions", xlab="RVSP_mmHg", cex.lab=2, cex.axis=1.5,
cex.main=2.5)
abline(a=fit.prop$coefficients[i,1], b=fit.prop$coefficients[i,2], col=4,
lwd=1)
}

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
rm(metadata_subset, props_subset, design, fit.prop)
design <- model.matrix(~ 0 + metadata$Group)
colnames(design) <- c("DCM", "HC")
mycontr <- makeContrasts(DCM-HC, levels=design)
differential_abundance_HC <- propeller.ttest(props, design, contrasts = mycontr, robust=TRUE, trend=FALSE, sort=TRUE)
differential_abundance_HC
PropMean.DCM PropMean.HC PropRatio Tstatistic P.Value FDR
CM3 0.11135354 0.03836732 2.9023016 3.800582 0.0006739595 0.003369797
CM4 0.08662428 0.02831785 3.0589987 3.202724 0.0032613900 0.008153475
CM1 0.50911696 0.57790022 0.8809773 -1.527054 0.1374288591 0.229048098
CM2 0.24032230 0.29095689 0.8259722 -1.002228 0.3244063073 0.405507884
CM5 0.05258291 0.06445772 0.8157738 -0.479165 0.6353640356 0.635364036
par(mfrow = c(1, 5), mar = c(3, 5, 3, 1))
for (i in seq(1, 5, 1)) {
stripchart(props$Proportions[i,] ~ metadata$Group,
vertical = TRUE, pch = 16, method = "jitter",
col = c("orange", "purple"),
cex = 2, ylab = "Proportions", cex.lab = 2, cex.axis = 2)
title(rownames(props$Proportions)[i], cex.main = 2.5)
}

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
rm(design, mycontr)
Remove genes expressed in <1% of cells and MT and RP genes
percent_stats <- Percent_Expressing(seurat_object = CM, features = rownames(CM), assay = "RNA", entire_object = TRUE)
percent_stats <- percent_stats %>%
dplyr::rename(percentCells = All_Cells)
percent_stats$gene <- rownames(percent_stats)
write.csv(percent_stats, file = here::here(output_dir_data, "CM_percent_stats.csv"), quote=F, row.names = T)
percent_stats <- read.csv(here::here(output_dir_data, "CM_percent_stats.csv"))
rownames(percent_stats) <- percent_stats$gene
keep_genes <- rownames(dplyr::filter(percent_stats, percentCells > 1))
data <- CM[which(rownames(CM) %in% keep_genes),]
rm(keep_genes)
t <- ezInteractiveTableRmd(percent_stats)
t[["x"]][["options"]][["pageLength"]] <- 10
t
Create and prepare pseudobulk object for differential expression analysis
pseudocounts <- Seurat2PB(data, sample="Sample", cluster = "cell_type")
colnames(pseudocounts) <- pseudocounts$samples$sample
saveRDS(pseudocounts,
here::here(output_dir_data, "CM_pseudocounts.rds"))
keep_samples <- pseudocounts$samples$lib.size > 5e4
pseudocounts <- pseudocounts[, keep_samples]
keep_genes <- filterByExpr(pseudocounts)
pseudocounts <- pseudocounts[keep_genes, ]
rm(keep_samples, keep_genes, data)
Create empty objects for the DE output
results <- list()
signif <- list()
volcano <- list()
n_de_genes <- data.frame()
pseudocounts_DCM <- pseudocounts[, metadata_DCM$Sample]
all.equal(metadata_DCM$Sample, colnames(pseudocounts_DCM$counts))
[1] TRUE
for (i in 6:(length(metadata_DCM))) {
metadata_subset <- metadata_DCM[complete.cases(metadata_DCM[[i]]),]
metadata_subset[,i] <- rescale(metadata_subset[,i])
design <- model.matrix(~ metadata_subset[[i]])
colnames(design)[2] <- colnames(metadata_subset[i])
pseudocounts_subset <- pseudocounts_DCM[, metadata_subset$Sample]
pseudocounts_subset <- normLibSizes(pseudocounts_subset)
pseudocounts_subset <- estimateDisp(pseudocounts_subset, design)
fit <- glmQLFit(pseudocounts_subset, design, robust=TRUE)
fit <- glmQLFTest(fit, coef = 2)
results[[i-5]] <- topTags(fit, n = Inf)$table
results[[i-5]] <- merge(results[[i-5]], percent_stats[, c("gene", "percentCells")], by = "gene", all.x = FALSE)
signif[[i-5]] <- results[[i-5]] %>% dplyr::filter(FDR < 0.05, abs(logFC) > 0.5) %>%
dplyr::arrange(FDR)
names(results)[i-5] <- colnames(metadata_subset[i])
names(signif)[i-5] <- colnames(metadata_subset[i])
write.csv(results[[i-5]], file = here::here(output_dir_data, paste0("CM_DE_Results_", colnames(metadata[i]), ".csv")), quote=F, row.names = F)
write.csv(signif[[i-5]], file = here::here(output_dir_data, paste0("CM_DE_Significant_", colnames(metadata[i]), ".csv")), quote=F, row.names = F)
n_de_genes_temp <- data.frame(metadata = colnames(metadata_subset[i]),
n_upregulated = sum(signif[[i-5]]$logFC > 0.5),
n_downregulated = sum(signif[[i-5]]$logFC < -0.5),
n_tot = nrow(signif[[i-5]])
)
n_de_genes <- rbind(n_de_genes, n_de_genes_temp)
volcano[[i-5]] <- EnhancedVolcano(results[[i-5]],
lab = results[[i-5]]$gene,
x = "logFC",
y = "FDR",
labSize = 0,
titleLabSize = 16,
subtitleLabSize = 14,
axisLabSize = 12,
captionLabSize = 9,
pointSize = 0.5,
FCcutoff = 0.5,
pCutoff = 0.05,
ylim = c(0, 3.5),
col = c("black", "black", "black", "red"),
colAlpha = 1,
drawConnectors = FALSE,
subtitle = NULL,
title = paste0("Cardiomyocytes", "\n", colnames(metadata_subset[i]), "\nn = ", nrow(metadata_subset))
) + theme(legend.position = "none")
names(volcano)[i-5] <- colnames(metadata_subset[i])
}
rm(metadata_subset, design, pseudocounts_subset, fit, n_de_genes_temp)
print((volcano[[1]] | volcano[[2]] | volcano[[3]] | volcano[[4]] | volcano[[5]] | volcano[[6]] | volcano[[7]]) /
(volcano[[8]] | volcano[[9]] | volcano[[10]] | volcano[[11]] | volcano[[12]] | volcano[[13]] | volcano[[14]]))

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
gene_sets <- list()
for (i in 1:length(signif)) {
gene_sets[[i]] <- signif[[i]]$gene
names(gene_sets)[i] <- names(signif)[i]
}
m <- make_comb_mat(gene_sets, mode = "intersect", min_set_size = 10)
m <- m[comb_degree(m) >= 2]
m <- m[comb_size(m) >= 10]
upset_plot <- UpSet(m, column_title = "Cardiomyocytes", comb_order = order(comb_size(m)),
top_annotation = upset_top_annotation(m, add_numbers = TRUE),
right_annotation = upset_right_annotation(m, add_numbers = TRUE))
draw(upset_plot)

# Calculate gene counts across cells
gene_counts <- rowSums(CM@assays$RNA@counts > 0)
# Filter genes to include only those expressed in at least 3 cells
universe_genes <- names(gene_counts[gene_counts >= 3])
write.csv(as.data.frame(universe_genes), here::here(output_dir_data, "CM_universe_genes.csv"), quote=F, row.names = F)
# Convert universe gene symbols to Entrez IDs
universe_entrez <- bitr(universe_genes, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)$ENTREZID
write.csv(as.data.frame(universe_entrez), here::here(output_dir_data, "CM_universe_entrez.csv"), quote=F, row.names = F)
rm(gene_counts)
up_genes <- filter(signif[["mPAP_mmHg"]], logFC > 0.5)$gene
up_genes <- bitr(up_genes, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")$ENTREZID
down_genes <- filter(signif[["mPAP_mmHg"]], logFC < -0.5)$gene
down_genes <- bitr(down_genes, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")$ENTREZID
reorder_GO_by_pvalue <- function(enrichGO_result) {
go_results <- enrichGO_result@result
go_results_sorted <- go_results[order(go_results$p.adjust), ]
enrichGO_result@result <- go_results_sorted
return(enrichGO_result)
}
mPAP_GO_up <- enrichGO(up_genes, OrgDb = "org.Hs.eg.db", universe = universe_entrez, ont = "ALL", readable = TRUE)
mPAP_GO_up <- gsfilter(mPAP_GO_up, by = 'Count', min = 3)
mPAP_GO_up <- reorder_GO_by_pvalue(mPAP_GO_up)
saveRDS(mPAP_GO_up, here::here(output_dir_data, "CM_ORA_GO_mPAP_up.rds"))
mPAP_GO_down <- enrichGO(down_genes, OrgDb = "org.Hs.eg.db", universe = universe_entrez, ont = "ALL", readable = TRUE)
mPAP_GO_down <- gsfilter(mPAP_GO_down, by = 'Count', min = 3)
mPAP_GO_down <- reorder_GO_by_pvalue(mPAP_GO_down)
saveRDS(mPAP_GO_down, here::here(output_dir_data, "CM_ORA_GO_mPAP_down.rds"))
mPAP_GO_up <- readRDS(here::here(output_dir_data, "CM_ORA_GO_mPAP_up.rds"))
mPAP_GO_down <- readRDS(here::here(output_dir_data, "CM_ORA_GO_mPAP_down.rds"))
p1 <- dotplot(mPAP_GO_up, showCategory = 10, title = paste0("GO - Over-representation analysis", "\nHigh mPAP - associated genes"), label_format = 27, font.size = 15) +
theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
axis.text.y = element_text(face = "bold"))
p2 <- dotplot(mPAP_GO_down, showCategory = 10, title = paste0("GO - Over-representation analysis", "\nLow mPAP - associated genes"), label_format = 27, font.size = 15) +
theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
axis.text.y = element_text(face = "bold"))
print(p1|p2)

mPAP_GO_simplified_up <- simplify(
mPAP_GO_up,
cutoff = 0.75,
by = "p.adjust",
select_fun = min,
measure = "Wang",
semData = NULL
)
mPAP_GO_simplified_down <- simplify(
mPAP_GO_down,
cutoff = 0.75,
by = "p.adjust",
select_fun = min,
measure = "Wang",
semData = NULL
)
p1 <- dotplot(mPAP_GO_simplified_up, showCategory = 10, title = paste0("GO - Over-representation analysis", "\nHigh mPAP - associated genes"), label_format = 27, font.size = 15) +
theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
axis.text.y = element_text(face = "bold"))
p2 <- dotplot(mPAP_GO_simplified_down, showCategory = 10, title = paste0("GO - Over-representation analysis", "\nLow mPAP - associated genes"), label_format = 27, font.size = 15) +
theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
axis.text.y = element_text(face = "bold"))
print(p1 | p2)

mPAP_KEGG_up <- enrichKEGG(up_genes, organism = 'hsa', universe = universe_entrez)
mPAP_KEGG_up <- gsfilter(mPAP_KEGG_up, by = 'Count', min = 3)
saveRDS(mPAP_KEGG_up, here::here(output_dir_data, "CM_ORA_KEGG_mPAP_up.rds"))
mPAP_KEGG_down <- enrichKEGG(down_genes, organism = 'hsa', universe = universe_entrez)
mPAP_KEGG_down <- gsfilter(mPAP_KEGG_down, by = 'Count', min = 3)
mPAP_KEGG_up <- readRDS(here::here(output_dir_data, "CM_ORA_KEGG_mPAP_up.rds"))
p1 <- dotplot(mPAP_KEGG_up, showCategory = 10, title = paste0("KEGG Pathways", "\nHigh mPAP - associated genes"), label_format = 27, font.size = 15) +
theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
axis.text.y = element_text(face = "bold"))
print(p1)

mPAP_REACTOME_up <- enrichPathway(up_genes, organism = 'human', universe = universe_entrez, readable = TRUE, pvalueCutoff = 0.1)
mPAP_REACTOME_up <- gsfilter(mPAP_REACTOME_up, by = 'Count', min = 3)
saveRDS(mPAP_REACTOME_up, here::here(output_dir_data, "CM_ORA_REACTOME_mPAP_up.rds"))
mPAP_REACTOME_down <- enrichPathway(down_genes, organism = 'human', universe = universe_entrez, readable = TRUE, pvalueCutoff = 0.1)
mPAP_REACTOME_down <- gsfilter(mPAP_REACTOME_down, by = 'Count', min = 3)
mPAP_REACTOME_up <- readRDS(here::here(output_dir_data, "CM_ORA_REACTOME_mPAP_up.rds"))
p1 <- dotplot(mPAP_REACTOME_up, showCategory = 10, title = paste0("REACTOME", "\nHigh mPAP - associated genes"), label_format = 27, font.size = 15) +
theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
axis.text.y = element_text(face = "bold"))
print(p1)

rm(mPAP_GO_up, mPAP_GO_down, mPAP_GO_simplified_up, mPAP_GO_simplified_down, mPAP_KEGG_up, mPAP_KEGG_down, mPAP_REACTOME_up, mPAP_REACTOME_down)
ranked_genes <- results[["mPAP_mmHg"]]
ranked_genes$metric <- ranked_genes$logFC*-log10(ranked_genes$PValue)
entrezid <- bitr(ranked_genes$gene, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db", drop = TRUE)
ranked_genes <- merge(ranked_genes, entrezid, by.x = "gene", by.y = "SYMBOL")
ranked_genes <- ranked_genes[!duplicated(ranked_genes$metric), ]
genelist <- ranked_genes$metric
names(genelist) <- ranked_genes$ENTREZID
genelist <- genelist[order(genelist, decreasing = T)]
write.csv(ranked_genes, file = here::here(output_dir_data, "CM_ranked_genes.csv"), quote=F, row.names = F)
t <- ezInteractiveTableRmd(ranked_genes)
t[["x"]][["options"]][["pageLength"]] <- 10
t
gseaGO_all <- gseGO(geneList = genelist,
OrgDb = org.Hs.eg.db,
ont = "ALL",
keyType = "ENTREZID",
nPermSimple = 1000000,
minGSSize = 10)
gseaGO_all <- reorder_GO_by_pvalue(gseaGO_all)
saveRDS(gseaGO_all, here::here(output_dir_data, "CM_GSEA_GO_mPAP.rds"))
gseaGO_all <- readRDS(here::here(output_dir_data, "CM_GSEA_GO_mPAP.rds"))
gseaGO_all_simplified <- simplify(
gseaGO_all,
cutoff = 0.75,
by = "p.adjust",
select_fun = min,
measure = "Wang",
semData = NULL
)
p1 <- dotplot(gseaGO_all_simplified, showCategory = 10, split=".sign", title = paste0("GO Enrichment Analysis", "\n mPAP-associated genes"), label_format = 30, font.size = 15) +
theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
axis.text.y = element_text(face = "bold"),
strip.text = element_text(face = "bold", size = 18)) +
facet_grid(.~.sign)
print(p1)

classify_terms <- function(df) {
df$sign <- ifelse(df$NES > 0, "activated", "suppressed")
return(df)
}
gseaGO_all_simplified@result <- classify_terms(gseaGO_all_simplified@result)
# Split the GSEA results
activated_terms <- subset(gseaGO_all_simplified@result, sign == "activated")
suppressed_terms <- subset(gseaGO_all_simplified@result, sign == "suppressed")
# Create two separate enrichResult objects
gseaGO_all_activated <- gseaGO_all_simplified
gseaGO_all_activated@result <- activated_terms
gseaGO_all_activated <- pairwise_termsim(gseaGO_all_activated)
gseaGO_all_suppressed <- gseaGO_all_simplified
gseaGO_all_suppressed@result <- suppressed_terms
gseaGO_all_suppressed <- pairwise_termsim(gseaGO_all_suppressed)
# Plot activated terms
p1 <- emapplot(gseaGO_all_activated, showCategory = 100, cex_label_category = 0.7, min_edge = 0.15) +
ggtitle("Activated Terms") +
theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5))
# Plot suppressed terms
p2 <- emapplot(gseaGO_all_suppressed, showCategory = 100, cex_label_category = 0.7, min_edge = 0.15) +
ggtitle("Suppressed Terms") +
theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5))
# Combine the plots
combined_plot <- p1 | p2
print(combined_plot)

p1 <- ridgeplot(gseaGO_all_simplified, showCategory = 20) + labs(x = "enrichment distribution")
print(p1)

gseaKEGG <- gseKEGG(geneList = genelist,
organism = "hsa",
keyType = "kegg",
nPermSimple = 100000)
saveRDS(gseaKEGG, here::here(output_dir_data, "CM_GSEA_KEGG_mPAP.rds"))
gseaKEGG <- readRDS(here::here(output_dir_data, "CM_GSEA_KEGG_mPAP.rds"))
p1 <- dotplot(gseaKEGG, showCategory = 10, split=".sign", title = paste0("GSEA\nKEGG Pathways", "\n mPAP-associated genes"), label_format = 30, font.size = 15) +
theme(plot.title = element_text( face = "bold", size = 18, hjust = 0.5),
axis.text.y = element_text(face = "bold"),
strip.text = element_text(face = "bold", size = 18)) +
facet_grid(.~.sign)
print(p1)

gseaKEGG@result <- classify_terms(gseaKEGG@result)
# Split the GSEA results
activated_terms <- subset(gseaKEGG@result, sign == "activated")
suppressed_terms <- subset(gseaKEGG@result, sign == "suppressed")
# Create two separate enrichResult objects
gseaKEGG_activated <- gseaKEGG
gseaKEGG_activated@result <- activated_terms
gseaKEGG_suppressed <- gseaKEGG
gseaKEGG_suppressed@result <- suppressed_terms
# Plot activated terms
p1 <- emapplot(pairwise_termsim(gseaKEGG_activated), showCategory = 20, cex_label_category = 0.7, min_edge = 0.15) +
ggtitle("Activated Terms") +
theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5))
# Plot suppressed terms
p2 <- emapplot(pairwise_termsim(gseaKEGG_suppressed), showCategory = 20, cex_label_category = 0.7, min_edge = 0.15) +
ggtitle("Suppressed Terms") +
theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5))
# Combine the plots
combined_plot <- p1 | p2
print(combined_plot)

p1 <- ridgeplot(gseaKEGG, showCategory = 20) + labs(x = "enrichment distribution")
print(p1)

gseaGO_all_activated_df <- setReadable(gseaGO_all_activated, OrgDb = org.Hs.eg.db)
gseaGO_all_activated_df <- gseaGO_all_activated_df@result
metadata_mPAP <- metadata[complete.cases(metadata[["mPAP_mmHg"]]),]
pseudocounts_mPAP <- pseudocounts[, metadata_mPAP$Sample]
pseudocounts_mPAP <- normLibSizes(pseudocounts_mPAP)
counts_mPAP <- cpm(pseudocounts_mPAP, log = TRUE)
## Set a color-blind friendly palette
heat_colors <- rev(brewer.pal(11, "PuOr"))
metadata_heatmap <- metadata_mPAP %>% dplyr::select(Sample, mPAP_mmHg, Batch_processing, Age_years) %>% dplyr::arrange(mPAP_mmHg)
DEGs_mPAP <- signif[["mPAP_mmHg"]]
sig_counts <- counts_mPAP[rownames(counts_mPAP) %in% DEGs_mPAP$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]
## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts,
color = heat_colors,
cluster_rows = TRUE,
show_rownames = FALSE,
annotation_col = metadata_heatmap,
border_color = NA,
fontsize = 10,
scale = "row",
cluster_cols=F,
fontsize_row = 10,
height = 20)

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
rm(metadata_heatmap, sig_counts)
metadata_heatmap <- metadata_mPAP %>% dplyr::select(Sample, mPAP_mmHg, Batch_processing, Age_years) %>% dplyr::arrange(mPAP_mmHg)
DEGs_mPAP_top <- DEGs_mPAP
DEGs_mPAP_top <- DEGs_mPAP_top %>%
dplyr::filter(percentCells > 5)
DEGs_mPAP_top_up <- DEGs_mPAP_top %>%
dplyr::filter(logFC > 2) %>%
dplyr::arrange(FDR)
DEGs_mPAP_top_down <- DEGs_mPAP_top %>%
dplyr::filter(logFC < -2) %>%
dplyr::arrange(FDR)
DEGs_mPAP_top <- dplyr::bind_rows(head(DEGs_mPAP_top_up, 20), head(DEGs_mPAP_top_down, 20))
rm(DEGs_mPAP_top_up, DEGs_mPAP_top_down)
sig_counts <- counts_mPAP[rownames(counts_mPAP) %in% DEGs_mPAP_top$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]
## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts,
color = heat_colors,
cluster_rows = TRUE,
annotation_col = metadata_heatmap,
border_color = NA,
fontsize = 10,
scale = "row",
cluster_cols=F,
fontsize_row = 15,
height = 20,
show_rownames = T)

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
rm(metadata_heatmap, sig_counts)
metadata_heatmap <- metadata %>%
dplyr::filter(Sample %in% metadata_mPAP$Sample | Group == "HC")
metadata_heatmap <- metadata_heatmap %>%
mutate(Group = factor(Group, levels = c("HC", "DCM"))) %>%
arrange(Group, mPAP_mmHg) %>%
dplyr::select(Sample, Group, mPAP_mmHg, Batch_processing, Age_years)
pseudocounts_mPAP_HC <- pseudocounts[, metadata_heatmap$Sample]
pseudocounts_mPAP_HC <- normLibSizes(pseudocounts_mPAP_HC)
counts_mPAP_HC <- cpm(pseudocounts_mPAP_HC, log = TRUE)
sig_counts <- counts_mPAP_HC[rownames(counts_mPAP_HC) %in% DEGs_mPAP$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]
## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts,
color = heat_colors,
cluster_rows = TRUE,
show_rownames = FALSE,
annotation_col = metadata_heatmap,
border_color = NA,
fontsize = 10,
scale = "row",
cluster_cols=F,
fontsize_row = 10,
height = 20)

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
rm(metadata_heatmap, sig_counts)
metadata_heatmap <- metadata %>%
dplyr::filter(Sample %in% metadata_mPAP$Sample | Group == "HC")
metadata_heatmap <- metadata_heatmap %>%
mutate(Group = factor(Group, levels = c("HC", "DCM"))) %>%
arrange(Group, mPAP_mmHg) %>%
dplyr::select(Sample, Group, mPAP_mmHg, Batch_processing, Age_years)
pseudocounts_mPAP_HC <- pseudocounts[, metadata_heatmap$Sample]
pseudocounts_mPAP_HC <- normLibSizes(pseudocounts_mPAP_HC)
counts_mPAP_HC <- cpm(pseudocounts_mPAP_HC, log = TRUE)
sig_counts <- counts_mPAP_HC[rownames(counts_mPAP_HC) %in% DEGs_mPAP_top$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]
## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts,
color = heat_colors,
cluster_rows = TRUE,
annotation_col = metadata_heatmap,
border_color = NA,
fontsize = 10,
scale = "row",
cluster_cols=F,
fontsize_row = 15,
height = 20,
show_rownames = T)

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
# Filter and arrange metadata
metadata_mPAP_HC <- metadata %>%
dplyr::filter(Sample %in% metadata_mPAP$Sample | Group == "HC") %>%
dplyr::mutate(Group = factor(Group, levels = c("HC", "DCM"))) %>%
dplyr::arrange(Group, mPAP_mmHg) %>%
dplyr::select(Sample, Group, mPAP_mmHg, Batch_processing, Age_years)
# Normalize pseudocounts and calculate log CPM
pseudocounts_mPAP_HC <- pseudocounts[, metadata_mPAP_HC$Sample]
pseudocounts_mPAP_HC <- normLibSizes(pseudocounts_mPAP_HC)
counts_mPAP_HC <- edgeR::cpm(pseudocounts_mPAP_HC, log = TRUE)
# Filter counts matrix by DEGs
counts_matrix <- counts_mPAP_HC[rownames(counts_mPAP_HC) %in% DEGs_mPAP$gene, ]
# Convert to long format and add Condition dividing DCM samples in 3 groups based on mPAP
counts_tibble <- counts_matrix %>%
as.data.frame() %>%
tibble::rownames_to_column(var = "gene") %>%
tidyr::pivot_longer(-gene, names_to = "sample", values_to = "log_pseudocounts") %>%
dplyr::mutate(Condition = dplyr::case_when(
sample %in% c("CZD4", "CZD5", "CZD6", "CZD7") ~ "HC",
sample %in% c("PK439", "PK375", "PK229", "P5", "P8", "PK153") ~ "DCM - Mild PH",
sample %in% c("PK385", "PK437", "PK323", "PK403", "PK213", "PK461") ~ "DCM - Moderate PH",
sample %in% c("PK369", "PK357", "PK337", "PK345", "PK401", "PK463") ~ "DCM - Severe PH"
))
# Compute the mean log counts for each gene and condition
mean_tibble <- counts_tibble %>%
dplyr::group_by(gene, Condition) %>%
dplyr::summarise(mean_log_counts = mean(log_pseudocounts, na.rm = TRUE), .groups = 'drop')
# Convert to wide format
wide_mean_tibble <- mean_tibble %>%
tidyr::pivot_wider(names_from = Condition, values_from = mean_log_counts)
# Convert the wide tibble to a matrix and transpose
mean_matrix <- wide_mean_tibble %>%
tibble::column_to_rownames(var = "gene") %>%
as.matrix()
mean_matrix <- mean_matrix %>%
t() %>%
scale() %>%
t()
# Perform hierarchical clustering
gene_dist <- dist(mean_matrix)
gene_hclust <- hclust(gene_dist, method = "complete")
# Plot the dendrogram
plot(gene_hclust, labels = FALSE)
abline(h = 3, col = "brown", lwd = 2) # Add horizontal line to illustrate cutting dendrogram

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
# Cut the dendrogram into 4 clusters
gene_cluster <- cutree(gene_hclust, k = 4) %>%
tibble::enframe() %>%
dplyr::rename(gene = name, cluster = value)
# Convert scaled matrix to long format and join with gene_cluster
mean_tibble_cluster <- mean_matrix %>%
as.data.frame() %>%
tibble::rownames_to_column(var = "gene") %>%
tidyr::pivot_longer(-gene, names_to = "Condition", values_to = "scaled_log_pseudocounts") %>%
dplyr::inner_join(gene_cluster, by = "gene")
# Set factor levels for Condition
mean_tibble_cluster$Condition <- factor(mean_tibble_cluster$Condition, levels = c("HC", "DCM - Mild PH", "DCM - Moderate PH", "DCM - Severe PH"))
# Calculate the number of genes per cluster
cluster_counts <- mean_tibble_cluster %>%
dplyr::group_by(cluster) %>%
dplyr::summarise(gene_count = dplyr::n_distinct(gene))
# Create a named vector for custom facet labels
custom_labels <- setNames(paste("Cluster", cluster_counts$cluster, "-", cluster_counts$gene_count, "genes"), cluster_counts$cluster)
# Plot the data with custom facet labels
mean_tibble_cluster %>%
ggplot(aes(Condition, scaled_log_pseudocounts)) +
geom_line(aes(group = gene), alpha = 0.2) +
geom_line(stat = "summary", fun = "median", colour = "brown", size = 1.5, aes(group = 1)) +
facet_grid(cols = vars(cluster), labeller = as_labeller(custom_labels)) +
theme(
strip.text = element_text(face = "bold", size = 12), # Make facet labels bold
axis.text.x = element_text(angle = 45, hjust = 1, size = 10) # Rotate x-axis labels by 45 degrees
)

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
DEGs_clusters <- list()
DEGs_clusters[["cluster1"]] <- unique(dplyr::filter(mean_tibble_cluster, cluster == "1")$gene)
DEGs_clusters[["cluster2"]] <- unique(dplyr::filter(mean_tibble_cluster, cluster == "2")$gene)
DEGs_clusters[["cluster3"]] <- unique(dplyr::filter(mean_tibble_cluster, cluster == "3")$gene)
DEGs_clusters[["cluster4"]] <- unique(dplyr::filter(mean_tibble_cluster, cluster == "4")$gene)
DEGs_clusters_entrez <- list()
DEGs_clusters_entrez[["cluster1"]] <- bitr(DEGs_clusters[["cluster1"]], fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")$ENTREZID
DEGs_clusters_entrez[["cluster2"]] <- bitr(DEGs_clusters[["cluster2"]], fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")$ENTREZID
DEGs_clusters_entrez[["cluster3"]] <- bitr(DEGs_clusters[["cluster3"]], fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")$ENTREZID
DEGs_clusters_entrez[["cluster4"]] <- bitr(DEGs_clusters[["cluster4"]], fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")$ENTREZID
cluster1_GO <- enrichGO(DEGs_clusters_entrez[["cluster1"]], OrgDb = "org.Hs.eg.db", universe = universe_entrez, ont = "ALL", readable = TRUE)
cluster1_GO <- gsfilter(cluster1_GO, by = 'Count', min = 3)
nrow(cluster1_GO)
cluster2_GO <- enrichGO(DEGs_clusters_entrez[["cluster2"]], OrgDb = "org.Hs.eg.db", universe = universe_entrez, ont = "ALL", readable = TRUE)
cluster2_GO <- gsfilter(cluster2_GO, by = 'Count', min = 3)
nrow(cluster2_GO)
cluster3_GO <- enrichGO(DEGs_clusters_entrez[["cluster3"]], OrgDb = "org.Hs.eg.db", universe = universe_entrez, ont = "ALL", readable = TRUE)
cluster3_GO <- gsfilter(cluster3_GO, by = 'Count', min = 10)
nrow(cluster3_GO)
cluster3_GO <- simplify(
cluster3_GO,
cutoff = 0.7,
by = "p.adjust",
select_fun = min,
measure = "Wang",
semData = NULL
)
cluster3_GO <- reorder_GO_by_pvalue(cluster3_GO)
saveRDS(cluster3_GO, here::here(output_dir_data, "CM_cluster3_ORA.rds"))
cluster3_GO <- readRDS(here::here(output_dir_data, "CM_cluster3_ORA.rds"))
dotplot(cluster3_GO,
showCategory = 15,
title = paste0("GO - Over-representation analysis", "\nDEGs cluster 3"),
label_format = 27,
font.size = 15) +
theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5),
axis.text.y = element_text(face = "bold"))

cluster4_GO <- enrichGO(DEGs_clusters_entrez[["cluster4"]], OrgDb = "org.Hs.eg.db", universe = universe_entrez, ont = "ALL", readable = TRUE)
cluster4_GO <- gsfilter(cluster4_GO, by = 'Count', min = 10)
nrow(cluster4_GO)
cluster4_GO <- simplify(
cluster4_GO,
cutoff = 0.7,
by = "p.adjust",
select_fun = min,
measure = "Wang",
semData = NULL
)
cluster4_GO <- reorder_GO_by_pvalue(cluster4_GO)
saveRDS(cluster4_GO, here::here(output_dir_data, "CM_cluster4_ORA.rds"))
cluster4_GO <- readRDS(here::here(output_dir_data, "CM_cluster4_ORA.rds"))
dotplot(cluster4_GO,
showCategory = 15,
title = paste0("GO - Over-representation analysis", "\nDEGs cluster 4"),
label_format = 27,
font.size = 15) +
theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5),
axis.text.y = element_text(face = "bold"))

cluster3_GO <- pairwise_termsim(cluster3_GO)
cluster4_GO <- pairwise_termsim(cluster4_GO)
p1 <- emapplot(cluster3_GO, showCategory = 100, cex_label_category = 0.7, min_edge = 0.15) +
ggtitle("DEGs cluster 3") +
theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5))
# Plot suppressed terms
p2 <- emapplot(cluster4_GO, showCategory = 100, cex_label_category = 0.7, min_edge = 0.15) +
ggtitle("DEGs cluster 4") +
theme(plot.title = element_text(face = "bold", size = 18, hjust = 0.5))
# Combine the plots
combined_plot <- p1 | p2
print(combined_plot)

DEGs_mPAP_cluster1 <- DEGs_mPAP %>%
filter(gene %in% DEGs_clusters[["cluster1"]])
top_DEGs_mPAP_cluster1 <- DEGs_mPAP_cluster1 %>%
filter(percentCells > 5) %>%
filter(FDR < 0.01) %>%
arrange(logFC) %>%
head(20)
DEGs_mPAP_cluster2 <- DEGs_mPAP %>%
filter(gene %in% DEGs_clusters[["cluster2"]])
top_DEGs_mPAP_cluster2 <- DEGs_mPAP_cluster2 %>%
filter(percentCells > 5) %>%
filter(FDR < 0.01) %>%
arrange(logFC) %>%
head(20)
DEGs_mPAP_cluster3 <- DEGs_mPAP %>%
filter(gene %in% DEGs_clusters[["cluster3"]])
top_DEGs_mPAP_cluster3 <- DEGs_mPAP_cluster3 %>%
filter(percentCells > 5) %>%
filter(FDR < 0.01) %>%
arrange(desc(logFC)) %>%
head(20)
DEGs_mPAP_cluster4 <- DEGs_mPAP %>%
filter(gene %in% DEGs_clusters[["cluster4"]])
top_DEGs_mPAP_cluster4 <- DEGs_mPAP_cluster4 %>%
filter(percentCells > 5) %>%
filter(FDR < 0.01) %>%
arrange(desc(logFC)) %>%
head(20)
metadata_heatmap <- metadata %>%
dplyr::filter(Sample %in% metadata_mPAP$Sample | Group == "HC")
metadata_heatmap <- metadata_heatmap %>%
mutate(Group = factor(Group, levels = c("HC", "DCM"))) %>%
arrange(Group, mPAP_mmHg) %>%
dplyr::select(Sample, Group, mPAP_mmHg, Age_years)
sig_counts <- counts_mPAP_HC[top_DEGs_mPAP_cluster1$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]
## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts,
color = heat_colors,
cluster_rows = F,
annotation_col = metadata_heatmap,
border_color = NA,
fontsize = 7,
scale = "row",
cluster_cols=F,
fontsize_row = 11,
cellheight = 17,
show_rownames = T,
main = "Top cluster 1 DEGs")

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
ezInteractiveTableRmd(top_DEGs_mPAP_cluster1)
metadata_heatmap <- metadata %>%
dplyr::filter(Sample %in% metadata_mPAP$Sample | Group == "HC")
metadata_heatmap <- metadata_heatmap %>%
mutate(Group = factor(Group, levels = c("HC", "DCM"))) %>%
arrange(Group, mPAP_mmHg) %>%
dplyr::select(Sample, Group, mPAP_mmHg, Age_years)
sig_counts <- counts_mPAP_HC[top_DEGs_mPAP_cluster2$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]
## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts,
color = heat_colors,
cluster_rows = F,
annotation_col = metadata_heatmap,
border_color = NA,
fontsize = 7,
scale = "row",
cluster_cols=F,
fontsize_row = 11,
cellheight = 17,
show_rownames = T,
main = "Top cluster 2 DEGs")

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
ezInteractiveTableRmd(top_DEGs_mPAP_cluster2)
metadata_heatmap <- metadata %>%
dplyr::filter(Sample %in% metadata_mPAP$Sample | Group == "HC")
metadata_heatmap <- metadata_heatmap %>%
mutate(Group = factor(Group, levels = c("HC", "DCM"))) %>%
arrange(Group, mPAP_mmHg) %>%
dplyr::select(Sample, Group, mPAP_mmHg, Age_years)
sig_counts <- counts_mPAP_HC[top_DEGs_mPAP_cluster3$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]
## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts,
color = heat_colors,
cluster_rows = F,
annotation_col = metadata_heatmap,
border_color = NA,
fontsize = 7,
scale = "row",
cluster_cols=F,
fontsize_row = 11,
cellheight = 17,
show_rownames = T,
main = "Top cluster 3 DEGs")

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
ezInteractiveTableRmd(top_DEGs_mPAP_cluster3)
metadata_heatmap <- metadata %>%
dplyr::filter(Sample %in% metadata_mPAP$Sample | Group == "HC")
metadata_heatmap <- metadata_heatmap %>%
mutate(Group = factor(Group, levels = c("HC", "DCM"))) %>%
arrange(Group, mPAP_mmHg) %>%
dplyr::select(Sample, Group, mPAP_mmHg, Age_years)
sig_counts <- counts_mPAP_HC[top_DEGs_mPAP_cluster4$gene, ]
sig_counts <- sig_counts[, metadata_heatmap$Sample, drop = FALSE]
## Run pheatmap using the metadata data frame for the annotation
pheatmap::pheatmap(sig_counts,
color = heat_colors,
cluster_rows = F,
annotation_col = metadata_heatmap,
border_color = NA,
fontsize = 7,
scale = "row",
cluster_cols=F,
fontsize_row = 11,
cellheight = 17,
show_rownames = T,
main = "Top cluster 4 DEGs")

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
ezInteractiveTableRmd(top_DEGs_mPAP_cluster4)
signif_cluster <- signif[["mPAP_mmHg"]] %>%
left_join(gene_cluster) %>%
left_join(ranked_genes)
write.csv(signif_cluster, file = here::here(output_dir_data, "CM_signif_cluster.csv"), quote=F, row.names = F)
ezInteractiveTableRmd(signif_cluster)
terms <- list(
muscle_structure_development = c("actin", "muscle", "sarcomere", "cardiac", "fibril", "fiber", "actomyosin"),
autophagy = c("autophag", "disassembly", "catabolic", "starvation"),
transcription = c("transcription", "receptor"),
vesicular_transport = c("transport", "Golgi", "vacuo", "vesicle", "endosom")
)
get_unique_genes <- function(df, terms) {
filtered_df <- df %>%
filter(grepl(paste(terms, collapse = "|"), Description, ignore.case = TRUE))
unique(unlist(strsplit(filtered_df$core_enrichment, split = "/")))
}
genes_of_interest <- lapply(terms, get_unique_genes, df = gseaGO_all_activated_df)
names(genes_of_interest) <- names(terms)
saveRDS(genes_of_interest, here::here(output_dir_data, "CM_genes_of_interest.rds"))
for (i in 1:length(genes_of_interest)) {
score_genes <- list(genes_of_interest[[i]])
score_name <- paste0(names(genes_of_interest)[i], "_score")
CM <- AddModuleScore(object = CM, features = score_genes, name = score_name)
}
metadata_CM <- CM@meta.data
colnames(metadata_CM) <- gsub("_score1", "_score", colnames(metadata_CM))
CM@meta.data <- metadata_CM
rm(metadata_CM)
plot_list <- lapply(paste0(names(genes_of_interest), "_score"), function(feature) {
p <- FeaturePlot(CM, features = feature, order = FALSE)
p + scale_color_viridis(discrete = FALSE, option = "turbo")
})
CombinePlots(plots = plot_list, ncol = 4)

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
plot_list <- lapply(paste0(names(genes_of_interest), "_score"), function(feature) {
p <- FeaturePlot(CM, features = feature, order = TRUE)
p + scale_color_viridis(discrete = FALSE, option = "turbo")
})
CombinePlots(plots = plot_list, ncol = 4)

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
VlnPlot(CM, features = paste0(names(genes_of_interest), "_score"), ncol = 4, group.by = "cell_state", pt.size = 0)

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
CM@meta.data[CM@meta.data$Sample %in% c("CZD4", "CZD5", "CZD6", "CZD7"), "Condition"] = "HC"
CM@meta.data[CM@meta.data$Sample %in% c("PK439", "PK375", "PK229", "P5", "P8", "PK153"), "Condition"] = "DCM - Mild PH"
CM@meta.data[CM@meta.data$Sample %in% c("PK385", "PK437", "PK323", "PK403", "PK213", "PK461"), "Condition"] = "DCM - Moderate PH"
CM@meta.data[CM@meta.data$Sample %in% c("PK369", "PK357", "PK337", "PK345", "PK401", "PK463"), "Condition"] = "DCM - Severe PH"
CM@meta.data[CM@meta.data$Sample %in% c("P3", "PK415", "PK441"), "Condition"] = "DCM - no mPAP data"
CM$Condition <- factor(CM$Condition, levels = c("HC", "DCM - Mild PH", "DCM - Moderate PH", "DCM - Severe PH", "DCM - no mPAP data"))
VlnPlot(CM, features = paste0(names(genes_of_interest), "_score"), ncol = 4, group.by = "Condition", pt.size = 0)

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
net <- get_collectri(organism='human', split_complexes=FALSE)
net
# A tibble: 43,178 × 3
source target mor
<chr> <chr> <dbl>
1 MYC TERT 1
2 SPI1 BGLAP 1
3 SMAD3 JUN 1
4 SMAD4 JUN 1
5 STAT5A IL2 1
6 STAT5B IL2 1
7 RELA FAS 1
8 WT1 NR0B1 1
9 NR0B2 CASP1 1
10 SP1 ALDOA 1
# ℹ 43,168 more rows
ranked_genes_TF <- ranked_genes %>%
dplyr::select(gene, logFC, metric, PValue) %>%
remove_rownames() %>%
column_to_rownames(var = "gene") %>%
as.matrix()
head(ranked_genes_TF)
logFC metric PValue
A1CF -4.4846666 -20.7823457 2.322259e-05
A2M -1.6943794 -6.0929917 2.535118e-04
A2M-AS1 0.2481086 0.1532492 2.411737e-01
A2ML1 -0.4535097 -0.1800881 4.007756e-01
A2ML1-AS1 -2.5607983 -8.7931153 3.683495e-04
A4GALT 2.0319539 6.1047409 9.899886e-04
# Run ulm
contrast_acts <- run_ulm(mat=ranked_genes_TF[, 'metric', drop=FALSE], net=net, .source='source', .target='target',
.mor='mor', minsize = 5)
head(contrast_acts)
# A tibble: 6 × 5
statistic source condition score p_value
<chr> <chr> <chr> <dbl> <dbl>
1 ulm ABL1 metric 1.49 0.137
2 ulm AHR metric 1.06 0.287
3 ulm AIP metric 0.893 0.372
4 ulm AIRE metric -0.730 0.465
5 ulm AP1 metric 0.776 0.438
6 ulm APEX1 metric 0.458 0.647
n_tfs <- 40
# Filter top TFs in both signs
f_contrast_acts <- contrast_acts %>%
mutate(rnk = NA)
msk <- f_contrast_acts$score > 0
f_contrast_acts[msk, 'rnk'] <- rank(-f_contrast_acts[msk, 'score'])
f_contrast_acts[!msk, 'rnk'] <- rank(-abs(f_contrast_acts[!msk, 'score']))
tfs <- f_contrast_acts %>%
arrange(rnk) %>%
head(n_tfs) %>%
pull(source)
f_contrast_acts <- f_contrast_acts %>%
filter(source %in% tfs)
# Plot
ggplot(f_contrast_acts, aes(x = reorder(source, score), y = score)) +
geom_bar(aes(fill = score), stat = "identity") +
scale_fill_gradient2(low = "darkblue", high = "indianred",
mid = "whitesmoke", midpoint = 0) +
theme_minimal() +
theme(axis.title = element_text(face = "bold", size = 12),
axis.text.x =
element_text(angle = 45, hjust = 1, size =10, face= "bold"),
axis.text.y = element_text(size =10, face= "bold"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
xlab("TFs")

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
FeaturePlot(CM, features = f_contrast_acts$source, ncol = 8)

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
contrast_acts_filtered <- contrast_acts %>%
filter(source %in% filter(percent_stats, percentCells > 5)$gene)
n_tfs <- 40
# Filter top TFs in both signs
f_contrast_acts <- contrast_acts_filtered %>%
mutate(rnk = NA) %>%
filter(p_value < 0.05)
msk <- f_contrast_acts$score > 0
f_contrast_acts[msk, 'rnk'] <- rank(-f_contrast_acts[msk, 'score'])
f_contrast_acts[!msk, 'rnk'] <- rank(-abs(f_contrast_acts[!msk, 'score']))
tfs <- f_contrast_acts %>%
arrange(rnk) %>%
head(n_tfs) %>%
pull(source)
f_contrast_acts <- f_contrast_acts %>%
filter(source %in% tfs)
# Plot
ggplot(f_contrast_acts, aes(x = reorder(source, score), y = score)) +
geom_bar(aes(fill = score), stat = "identity") +
scale_fill_gradient2(low = "darkblue", high = "indianred",
mid = "whitesmoke", midpoint = 0) +
theme_minimal() +
theme(axis.title = element_text(face = "bold", size = 12),
axis.text.x =
element_text(angle = 45, hjust = 1, size =10, face= "bold"),
axis.text.y = element_text(size =10, face= "bold"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
xlab("TFs")

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
FeaturePlot(CM, features = f_contrast_acts$source, ncol = 6)

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
tf <- 'MYOCD'
df <- net %>%
filter(source == tf) %>%
arrange(target) %>%
mutate(ID = target, color = "3") %>%
column_to_rownames('target')
inter <- sort(intersect(rownames(ranked_genes_TF),rownames(df)))
df <- df[inter, ]
df[,c('logfc', 't_value', 'p_value')] <- ranked_genes_TF[inter, ]
df <- df %>%
mutate(color = if_else(mor > 0 & t_value > 0, '1', color)) %>%
mutate(color = if_else(mor > 0 & t_value < 0, '2', color)) %>%
mutate(color = if_else(mor < 0 & t_value > 0, '2', color)) %>%
mutate(color = if_else(mor < 0 & t_value < 0, '1', color))
ggplot(df, aes(x = logfc, y = -log10(p_value), color = color, size=abs(mor))) +
geom_point() +
scale_colour_manual(values = c("red","royalblue3","grey")) +
geom_label_repel(aes(label = ID, size=1)) +
theme_minimal() +
theme(legend.position = "none") +
geom_vline(xintercept = 0, linetype = 'dotted') +
geom_hline(yintercept = 0, linetype = 'dotted') +
ggtitle(tf)

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
tf <- 'SRF'
df <- net %>%
filter(source == tf) %>%
arrange(target) %>%
mutate(ID = target, color = "3") %>%
column_to_rownames('target')
inter <- sort(intersect(rownames(ranked_genes_TF),rownames(df)))
df <- df[inter, ]
df[,c('logfc', 't_value', 'p_value')] <- ranked_genes_TF[inter, ]
df <- df %>%
mutate(color = if_else(mor > 0 & t_value > 0, '1', color)) %>%
mutate(color = if_else(mor > 0 & t_value < 0, '2', color)) %>%
mutate(color = if_else(mor < 0 & t_value > 0, '2', color)) %>%
mutate(color = if_else(mor < 0 & t_value < 0, '1', color))
ggplot(df, aes(x = logfc, y = -log10(p_value), color = color, size=abs(mor))) +
geom_point() +
scale_colour_manual(values = c("red","royalblue3","grey")) +
geom_label_repel(aes(label = ID, size=1)) +
theme_minimal() +
theme(legend.position = "none") +
geom_vline(xintercept = 0, linetype = 'dotted') +
geom_hline(yintercept = 0, linetype = 'dotted') +
ggtitle(tf)

tf <- 'MEF2A'
df <- net %>%
filter(source == tf) %>%
arrange(target) %>%
mutate(ID = target, color = "3") %>%
column_to_rownames('target')
inter <- sort(intersect(rownames(ranked_genes_TF),rownames(df)))
df <- df[inter, ]
df[,c('logfc', 't_value', 'p_value')] <- ranked_genes_TF[inter, ]
df <- df %>%
mutate(color = if_else(mor > 0 & t_value > 0, '1', color)) %>%
mutate(color = if_else(mor > 0 & t_value < 0, '2', color)) %>%
mutate(color = if_else(mor < 0 & t_value > 0, '2', color)) %>%
mutate(color = if_else(mor < 0 & t_value < 0, '1', color))
ggplot(df, aes(x = logfc, y = -log10(p_value), color = color, size=abs(mor))) +
geom_point() +
scale_colour_manual(values = c("red","royalblue3","grey")) +
geom_label_repel(aes(label = ID, size=1)) +
theme_minimal() +
theme(legend.position = "none") +
geom_vline(xintercept = 0, linetype = 'dotted') +
geom_hline(yintercept = 0, linetype = 'dotted') +
ggtitle(tf)

https://www.nature.com/articles/s41419-023-05665-8
tf <- 'MEF2C'
df <- net %>%
filter(source == tf) %>%
arrange(target) %>%
mutate(ID = target, color = "3") %>%
column_to_rownames('target')
inter <- sort(intersect(rownames(ranked_genes_TF),rownames(df)))
df <- df[inter, ]
df[,c('logfc', 't_value', 'p_value')] <- ranked_genes_TF[inter, ]
df <- df %>%
mutate(color = if_else(mor > 0 & t_value > 0, '1', color)) %>%
mutate(color = if_else(mor > 0 & t_value < 0, '2', color)) %>%
mutate(color = if_else(mor < 0 & t_value > 0, '2', color)) %>%
mutate(color = if_else(mor < 0 & t_value < 0, '1', color))
ggplot(df, aes(x = logfc, y = -log10(p_value), color = color, size=abs(mor))) +
geom_point() +
scale_colour_manual(values = c("red","royalblue3","grey")) +
geom_label_repel(aes(label = ID, size=1)) +
theme_minimal() +
theme(legend.position = "none") +
geom_vline(xintercept = 0, linetype = 'dotted') +
geom_hline(yintercept = 0, linetype = 'dotted') +
ggtitle(tf)

tf <- 'NFATC3'
df <- net %>%
filter(source == tf) %>%
arrange(target) %>%
mutate(ID = target, color = "3") %>%
column_to_rownames('target')
inter <- sort(intersect(rownames(ranked_genes_TF),rownames(df)))
df <- df[inter, ]
df[,c('logfc', 't_value', 'p_value')] <- ranked_genes_TF[inter, ]
df <- df %>%
mutate(color = if_else(mor > 0 & t_value > 0, '1', color)) %>%
mutate(color = if_else(mor > 0 & t_value < 0, '2', color)) %>%
mutate(color = if_else(mor < 0 & t_value > 0, '2', color)) %>%
mutate(color = if_else(mor < 0 & t_value < 0, '1', color))
ggplot(df, aes(x = logfc, y = -log10(p_value), color = color, size=abs(mor))) +
geom_point() +
scale_colour_manual(values = c("red","royalblue3","grey")) +
geom_label_repel(aes(label = ID, size=1)) +
theme_minimal() +
theme(legend.position = "none") +
geom_vline(xintercept = 0, linetype = 'dotted') +
geom_hline(yintercept = 0, linetype = 'dotted') +
ggtitle(tf)

# Run ulm
sample_acts <- run_ulm(mat=counts_mPAP_HC, net=net, .source='source', .target='target',
.mor='mor', minsize = 5)
# Transform to wide matrix
sample_acts_mat <- sample_acts %>%
pivot_wider(id_cols = 'condition', names_from = 'source',
values_from = 'score') %>%
column_to_rownames('condition') %>%
as.matrix()
# Get top tfs with more variable means across clusters
tfs <- arrange(f_contrast_acts, desc(score))$source
sample_acts_mat <- sample_acts_mat[,tfs]
# Scale per sample
sample_acts_mat <- scale(sample_acts_mat)
sample_acts_mat <- t(sample_acts_mat)
sample_acts_mat <- sample_acts_mat[, metadata_mPAP_HC$Sample, drop = FALSE]
# Choose color palette
palette_length = 100
my_color = colorRampPalette(c("Darkblue", "white","red"))(palette_length)
my_breaks <- c(seq(-3, 0, length.out=ceiling(palette_length/2) + 1),
seq(0.05, 3, length.out=floor(palette_length/2)))
# Plot
pheatmap::pheatmap(sample_acts_mat,
border_color = NA,
color=my_color,
breaks = my_breaks,
cluster_rows = F,
cluster_cols = F,
annotation_col = dplyr::select(metadata_mPAP_HC, -Batch_processing))

| Version | Author | Date |
|---|---|---|
| 61e8d17 | GinoBonazza | 2024-08-01 |
saveRDS(CM,
here::here(output_dir_data, "CM.rds"))
sessionInfo()
R version 4.3.1 (2023-06-16)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.3 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=en_US.UTF-8
[9] LC_ADDRESS=en_US.UTF-8 LC_TELEPHONE=en_US.UTF-8
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=en_US.UTF-8
time zone: Etc/UTC
tzcode source: system (glibc)
attached base packages:
[1] grid stats4 stats graphics grDevices utils datasets
[8] methods base
other attached packages:
[1] htmltools_0.5.8.1 DT_0.32
[3] ComplexHeatmap_2.18.0 UpSetR_1.4.0
[5] ezRun_3.19.1 Biostrings_2.70.2
[7] XVector_0.42.0 enrichplot_1.22.0
[9] dorothea_1.14.1 OmnipathR_3.10.1
[11] decoupleR_2.9.7 EnhancedVolcano_1.20.0
[13] ggrepel_0.9.5 speckle_1.2.0
[15] edgeR_4.0.16 limma_3.58.1
[17] ReactomePA_1.46.0 AnnotationHub_3.10.0
[19] BiocFileCache_2.10.1 dbplyr_2.4.0
[21] org.Hs.eg.db_3.18.0 AnnotationDbi_1.64.1
[23] clusterProfiler_4.10.1 gprofiler2_0.2.3
[25] harmony_1.2.0 clustree_0.5.1
[27] ggraph_2.2.1 rhdf5_2.46.1
[29] Rcpp_1.0.12 data.table_1.15.2
[31] plyr_1.8.9 gtable_0.3.4
[33] gtools_3.9.5 ArchR_1.0.2
[35] statmod_1.5.0 patchwork_1.2.0
[37] scater_1.30.1 scuttle_1.12.0
[39] Matrix.utils_0.9.7 RColorBrewer_1.1-3
[41] scales_1.3.0 knitr_1.45
[43] gridExtra_2.3 png_0.1-8
[45] pheatmap_1.0.12 reshape2_1.4.4
[47] lubridate_1.9.3 forcats_1.0.0
[49] stringr_1.5.1 purrr_1.0.2
[51] tidyr_1.3.1 tibble_3.2.1
[53] tidyverse_2.0.0 magrittr_2.0.3
[55] ggplot2_3.5.0 dplyr_1.1.4
[57] scCustomize_2.1.2 scDblFinder_1.16.0
[59] Matrix_1.6-5 DropletUtils_1.22.0
[61] SingleCellExperiment_1.24.0 SummarizedExperiment_1.32.0
[63] Biobase_2.62.0 GenomicRanges_1.54.1
[65] GenomeInfoDb_1.38.7 IRanges_2.36.0
[67] S4Vectors_0.40.2 BiocGenerics_0.48.1
[69] MatrixGenerics_1.14.0 matrixStats_1.2.0
[71] SeuratObject_5.0.2 Seurat_4.4.0
[73] xlsx_0.6.5 readxl_1.4.3
[75] readr_2.1.5 here_1.0.1
loaded via a namespace (and not attached):
[1] igraph_2.0.3 graph_1.80.0
[3] ica_1.0-3 plotly_4.10.4
[5] rematch2_2.1.2 zlibbioc_1.48.0
[7] tidyselect_1.2.1 rvest_1.0.4
[9] bit_4.0.5 doParallel_1.0.17
[11] clue_0.3-65 lattice_0.22-5
[13] rjson_0.2.21 blob_1.2.4
[15] S4Arrays_1.2.1 parallel_4.3.1
[17] cli_3.6.2 ggplotify_0.1.2
[19] goftest_1.2-3 BiocIO_1.12.0
[21] bluster_1.12.0 grr_0.9.5
[23] BiocNeighbors_1.20.2 uwot_0.1.16
[25] shadowtext_0.1.3 curl_5.2.1
[27] mime_0.12 evaluate_0.23
[29] tidytree_0.4.6 leiden_0.4.3.1
[31] stringi_1.8.3 backports_1.4.1
[33] XML_3.99-0.16.1 httpuv_1.6.14
[35] paletteer_1.6.0 rappdirs_0.3.3
[37] splines_4.3.1 logger_0.3.0
[39] bcellViper_1.38.0 sctransform_0.4.1
[41] ggbeeswarm_0.7.2 DBI_1.2.2
[43] HDF5Array_1.30.1 jquerylib_0.1.4
[45] reactome.db_1.86.2 withr_3.0.0
[47] git2r_0.33.0 rprojroot_2.0.4
[49] xgboost_1.7.7.1 lmtest_0.9-40
[51] ggnewscale_0.4.10 tidygraph_1.3.1
[53] rtracklayer_1.62.0 BiocManager_1.30.22
[55] htmlwidgets_1.6.4 fs_1.6.3
[57] labeling_0.4.3 SparseArray_1.2.4
[59] cellranger_1.1.0 reticulate_1.35.0
[61] zoo_1.8-12 timechange_0.3.0
[63] foreach_1.5.2 fansi_1.0.6
[65] ggtree_3.10.1 R.oo_1.26.0
[67] irlba_2.3.5.1 ggrastr_1.0.2
[69] gridGraphics_0.5-1 ellipsis_0.3.2
[71] lazyeval_0.2.2 yaml_2.3.8
[73] survival_3.5-8 scattermore_1.2
[75] BiocVersion_3.18.1 crayon_1.5.2
[77] RcppAnnoy_0.0.22 progressr_0.14.0
[79] tweenr_2.0.3 later_1.3.2
[81] ggridges_0.5.6 codetools_0.2-19
[83] GlobalOptions_0.1.2 KEGGREST_1.42.0
[85] Rtsne_0.17 shape_1.4.6.1
[87] Rsamtools_2.18.0 filelock_1.0.3
[89] pkgconfig_2.0.3 xml2_1.3.6
[91] GenomicAlignments_1.38.2 aplot_0.2.2
[93] spatstat.sparse_3.0-3 ape_5.7-1
[95] viridisLite_0.4.2 xtable_1.8-4
[97] highr_0.10 httr_1.4.7
[99] tools_4.3.1 globals_0.16.3
[101] beeswarm_0.4.0 checkmate_2.3.1
[103] nlme_3.1-164 selectr_0.4-2
[105] HDO.db_0.99.1 crosstalk_1.2.1
[107] digest_0.6.35 farver_2.1.1
[109] tzdb_0.4.0 yulab.utils_0.1.4
[111] viridis_0.6.5 glue_1.7.0
[113] cachem_1.0.8 polyclip_1.10-6
[115] generics_0.1.3 parallelly_1.37.1
[117] ScaledMatrix_1.10.0 pbapply_1.7-2
[119] vroom_1.6.5 spam_2.10-0
[121] gson_0.1.0 dqrng_0.3.2
[123] utf8_1.2.4 graphlayouts_1.1.1
[125] shiny_1.8.0 GenomeInfoDbData_1.2.11
[127] R.utils_2.12.3 rhdf5filters_1.14.1
[129] RCurl_1.98-1.14 memoise_2.0.1
[131] rmarkdown_2.26 R.methodsS3_1.8.2
[133] future_1.33.1 RANN_2.6.1
[135] Cairo_1.6-2 spatstat.data_3.0-4
[137] rstudioapi_0.15.0 cluster_2.1.6
[139] whisker_0.4.1 janitor_2.2.0
[141] spatstat.utils_3.0-4 hms_1.1.3
[143] fitdistrplus_1.1-11 munsell_0.5.0
[145] cowplot_1.1.3 colorspace_2.1-0
[147] rlang_1.1.4 DelayedMatrixStats_1.24.0
[149] sparseMatrixStats_1.14.0 dotCall64_1.1-1
[151] ggforce_0.4.2 circlize_0.4.16
[153] xfun_0.42 iterators_1.0.14
[155] abind_1.4-5 GOSemSim_2.28.1
[157] interactiveDisplayBase_1.40.0 treeio_1.26.0
[159] rJava_1.0-11 Rhdf5lib_1.24.2
[161] bitops_1.0-7 promises_1.2.1
[163] scatterpie_0.2.1 RSQLite_2.3.5
[165] qvalue_2.34.0 fgsea_1.28.0
[167] DelayedArray_0.28.0 GO.db_3.18.0
[169] compiler_4.3.1 prettyunits_1.2.0
[171] beachmat_2.18.1 graphite_1.48.0
[173] listenv_0.9.1 workflowr_1.7.1
[175] BiocSingular_1.18.0 tensor_1.5
[177] MASS_7.3-60.0.1 progress_1.2.3
[179] BiocParallel_1.36.0 spatstat.random_3.2-3
[181] R6_2.5.1 fastmap_1.1.1
[183] fastmatch_1.1-4 vipor_0.4.7
[185] ROCR_1.0-11 rsvd_1.0.5
[187] KernSmooth_2.23-22 miniUI_0.1.1.1
[189] deldir_2.0-4 bit64_4.0.5
[191] spatstat.explore_3.2-6 lifecycle_1.0.4
[193] ggprism_1.0.4 restfulr_0.0.15
[195] xlsxjars_0.6.1 sass_0.4.9
[197] vctrs_0.6.5 spatstat.geom_3.2-9
[199] snakecase_0.11.1 DOSE_3.28.2
[201] scran_1.30.2 ggfun_0.1.4
[203] sp_2.1-3 future.apply_1.11.1
[205] bslib_0.6.1 pillar_1.9.0
[207] metapod_1.10.1 locfit_1.5-9.9
[209] jsonlite_1.8.8 GetoptLong_1.0.5